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1  THE  BGK  MODEL 


1  The  BGK  model 


The  correct  description  of  rarefaction  effects  requires  replacing  hydrodynamic  equations 
with  the  more  general  Boltzmann  equation  (Cercignani,  1988)  which  takes  the  form: 

%  +  v-Vrf+-V„-(Ff)  =  C(f,f)  (1) 

at  m 

when  written  for  a  simple  monatomic  gas  composed  of  atoms  of  mass  m.  In  Eq.  (1),  the 
distribution  function  f(r,  v\t)  gives  the  number  of  atoms  with  position  r  and  velocity  v  at 
time  t.  The  vector  field  F{r,v\t)  describes  the  effects  of  an  assigned  external  force  field. 
The  rate  of  change  of  /  due  to  atomic  interactions  is  given  by  the  source  term 
which  is  a  non-linear  functional  of  /,  whose  precise  structure  depends  on  the  assumed 
atomic  interaction  forces.  For  instance,  hard  sphere  interaction  leads  to  the  following 
form  of  the  collision  integral: 

C(f,f)  =  %r  [  dvi  [  d2k  [f  (x ,  v\\t)  f  (x ,  v* \t)  -  /(*,«i|i)/(*,v|i)]  \k  ■  vr\  (2) 
1  Jb?  Js2 

In  Eq.  (2),  a  is  the  hard  sphere  diameter,  vr  =  V\  —  v  is  the  relative  velocity  of  two 
colliding  atoms  whereas  k  is  a  vector,  belonging  to  the  unit  sphere  S2,  used  to  specify  the 
relative  position  of  two  atoms  at  the  time  of  their  impact.  The  collisional  dynamics  of 
two  colliding  atoms  determines  the  pre-collisional  velocities  v*  and  v\  which  are  changed 
into  v  and  V\  by  a  binary  collision.  In  this  case  v*  and  v\  are  obtained  from  v ,  V\  and 
the  impact  vector  k  by  the  simple  relationships: 

v*  =  v  +  (vr  ■  k)k  v\  =  V\  —  (vr  ■  k)k  (3) 


The  complicated  structure  of  the  Boltzmann  collision  integral  leads  quite  naturally  to  ask 
weather  it  is  possible  to  replace  C(f,  f )  with  a  collision  term  which  has  a  simpler  mathe¬ 
matical  form  and  yet  preserves  the  basic  properties  of  the  original  expression.  The  simplest 
and  most  widely  used  kinetic  model  equation  has  been  proposed  in  1954  by  Bhatna- 
gar,  Gross  and  Krook(Bhatnagar  et  ah,  1954)  and  independently  by  Welander(Welander, 
1954).  In  the  BGK  model  C(f,f)  is  replaced  by  a  relaxation  term  having  the  following 
form: 

=  v[$(r,v\t)~  f(r,v\t)]  (4) 

BGKW 


In  Eq.(4),  v  is  the  collision  frequency,  whereas  $(r,u|t)  is  a  Maxwellian  distribution 
function: 


<E>(r,  v\t) 


n{r\t) 


2nRT(r\t ) 


3/2 


exp 


[  i  [v  - 

1  2  2RT(r\t)  J 


(5) 


characterized  by  the  density  n(r\t),  bulk  velocity  u[r\t)  and  temperature  T(r\t);  the 
gas  constant  R  is  defined  as  the  ratio  — ,  being  ks  the  Boltzmann  constant  and  m  the 
molecular  mass. 

The  form  of  the  collision  term  given  by  Eq.  (4)  can  somehow  be  justified  by  observing 
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1  THE  BGK  MODEL 


1.1  Basic  Model  Properties 


that  C(f,  f)  can  be  rewritten  as  the  difference  of  two  terms: 


C(fJ) 

"(r,v\t) 


G(f,  f )  -  v(r,v\t)f(r,v\t) 


TT  /  dvl 


drk 


I R3 


Is2 


f(x,vl\t)f{x,v*\t)\k  ■  vr 


a 

y 


<  R?  Js 2 


f(x,  Vi\t)\k  ■  vr\  dv !  d2fc 


(6) 

(7) 

(8) 


The  loss  term  z/(r,  v\t)f(r,  v\t )  has  exactly  the  same  form  both  in  the  original  and  in  the 
BGKW  model.  The  gain  term  in  the  Boltzmann  equation  describes  the  effects 

of  the  ’’arrival”  of  new  molecules  at  the  phase  space  point  (r,  v)  as  a  result  of  collisions. 
Since  collisions  tend  to  drive  the  gas  toward  a  local  equilibrium  condition,  described  by  a 
Maxwellian  distribution  function,  the  BGKW  model  assumes  that  collisions  immediately 
therm alize  molecules  through  the  approximate  gain  term  z/$(r,  v\ t).  The  model  structure 
is  made  much  simpler  by  assuming  that  the  collision  frequency  v  does  not  depend  on 
velocity  v ,  although  it  is  allowed  to  let  it  depend  on  local  macroscopic  quantities. 

If  collision  frequency  does  not  depend  on  velocity,  then  the  requirement  of  local  mass, 
momentum  and  energy  conservation  leads  to  the  relationships: 


I  &(v)dv  =  j  f(v)dv 
(  v$(v)  dv  —  J  vf{v )  dv 
v2${v)dv  =  [  v2f{v)dv 


which  determine  the  moments  n,  u  and  T  as: 

n  =  n  —  J  f{v)dv 

1  f 

u  =  u  —  —  vf(v )  dv 

n  J 


f  =  T  =  L  f  (v  -  uff(v)  dv 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 


In  other  words,  the  Maxwellian  density,  bulk  velocity  and  temperature  coincide  with  the 
corresponding  moments  of  the  distribution  function  /. 


1.1  Basic  Model  Properties 

In  spite  of  the  linear  apparence  of  the  BGKW  collision  term,  its  non-linearity  is  much 
stronger  (exponential)  than  the  quadratic  Boltzmann  collision  integral.  However,  the 
derivation  of  many  important  model  properties  turns  out  to  be  much  simpler. 

The  gas  homogeneous  relaxation  toward  the  final  equilibrium  state  is  described  by  the 
equation: 

^  =  v{n,  T)  [$(r,  v\t)  -  f(r,v\t)}  (17) 
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1  THE  BGK  MODEL 


Mass,  momentum  and  energy  conservations  keep  the  moments  n,  u  and  T,  constant 
during  the  evolution  of  /.  Hence,  their  constant  values  no,  u0  and  T0  are  determined  by 
the  initial  gas  state  /o  =  /(u|0).  Being  completely  determined  by  no,  u0  and  T0,  both 
the  Maxwellian  and  the  collision  frequency  n,  will  be  time  independent,  too.  Hence, 
the  solution  of  Eq.  (17)  can  easily  written  as  a  linear  combination  of  the  initial  and  final 
state  at  each  time  t: 

f(v\t)  =  f0(v)e~uot  +  (1  -  e~Uot)^(v)  (18) 


The  possibility  of  determining  a  simple  explicit  expression  for  the  homogeneous  relaxation 
is  quite  useful  in  the  application  of  fractional  time  step  schemes  to  the  numerical  solution 
of  BGKW  model  equation. 

In  the  space  homogeneous  case,  the  BGKW  model  also  allows  a  simple  proof  of  the  77- 
theorem  (Cercignani  ,  1988).  As  in  the  case  of  the  full  Boltzmann  equation,  the  entropy 
functional  is  dehned  as  'H(t)  =  f  f  (v)  log[f  (v\t)]  dv .  The  following  equality  is  easily 
obtained: 


^  =  1/  /(*  -  /)  log(/)  dv 


(19) 


Since  log(<f>)  is  a  linear  combination  of  the  collision  invariants,  the  r.h.s.  of  Eq.  (19)  can 
be  also  written  as 


=  v  /(*  -  /)  log(//4)  dv  +  v  /($  -  /)  log(4)  dv 


(20) 


being  the  last  integral  equal  to  zero.  The  remaining  term  gives 


(M 

dt 


<  0 


since  (<L  —  /)  log(//<L)  <  0,  the  equality  applying  only  if  /  =  $. 


(21) 


1.2  Hydrodynamic  limit  of  the  BGKW  Model 

Investigating  the  hydrodynamic  limit  of  the  BGKW  model  is  extremely  important  to 
judge  about  the  model  properties  and  capabilities  as  well  as  to  determine  the  proper  form 
of  the  collision  collision  frequency  v. 

The  classical  Chapman-Enskog  (Ferziger  and  Kaper,  1972)  expansion  provides  the  sim¬ 
plest  way  to  establish  a  connection  between  nu  and  the  transport  properties  predicted  by 
Eq.  (4).  It  is  convenient  to  rescale  the  space  variable  r  to  the  flow  macroscopic  reference 
scale  L  and  time  to  the  slow  time  scale  L/ dj RTref,  being  RTref  a  reference  thermal 
speed  value.  Since  v  is  of  the  order  of  — A ,  being  rref  a  reference  value  of  the  mean  free 
time,  the  rescaled  Eq.  (4)  reads 

e  +  v  ■  Vrf^j  =  v  [$(r,  v\t)  -  f(r,  v\ t)]  (22) 

where  the  (assumed)  small  parameter  e  =  \/Kr,,^T,'_L  pas  the  meaning  of  a  reference 
Knudsen  number.  In  the  limit  of  small  e,  collisions  dominate  the  flow  evolutions;  it  is 
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1.2  Hydrodynamic  limit  of  the  BGKW  Model 


then  assumed  that  kinetic  (i.e.  non  hydrodynamic  modes)  rapidly  decay  and  that  the 
distribution  function  is  determined  by  the  macroscopic  quantities 


' P(r\t) ' 

(3  =  |  u(r\t)  I  = 

KT(r\t)J  \^r 


m 

rnv 


m(v  —  u) 


f(r,v\t)dv 


(23) 


and  their  spatial  gradients.  Accordingly,  it  is  assumed  that 

f(v\(3,Vr(3,VrVr(3,...)  (24) 

The  macroscopic  fields  p,  u ,  and  T  appear  in  the  mass,  momentum  and  kinetic  energy 
conservation  equations: 

Dp 


P 


Cvp 


Dt 

Du 

~Dt 

DT 

dt 


— pV  o  u 

-Vo  P 

—V  o  q  +  P  :  Vv 


(25) 

(26) 
(27) 


where  Cv  =  M  is  the  gas  specific  heat  per  unit  mass,  whereas  the  kinetic  stress  tensor  P 
and  the  kinetic  heat  flux  vector  q  are  defined  as 


m(v  -  u)(v  -  u))\  ,  ,  . 

f(r,v\t)dv 


^{v  —  u)  [v  —  u 


(28) 


Taking  into  account  Eq.  (24),  the  conservation  equations  written  above  can  be  formally 
written  as 

33 

^  =  $(/ 3,VP/3,VPVr/3,...)  (29) 

The  dependence  of  /  on  macroscopic  helds  can  be  determined  from  the  expansion 

/  =  /(0)  +  e/(1)  +  e2/(2)  +  •  •  •  (30) 

which  also  induces  a  similar  expansion  of  the  streaming  term.  Hence,  Eq.  (22)  can  be 
formally  written  as 


D£ 

Dt 


(0) 


+  e 


Df\ 

Dt ) 


(i) 


+  . . . 


=  v 


($  -  /(«>  -  £/<‘>  -  ey<2)  + . . . ) 


(31) 


The  zero  and  first  order  terms  are  then  readily  obtained  by  equating  terms  of  correspond¬ 
ing  order  in  e: 


(0) 


/(0)  =  Hr,v\t) 

>"  --im 

where  the  term  (7^)^  is  obtained  from  the  expression 


(32) 

(33) 


<9/(°) 

~dT 


V-Vrf 


(0) 


(34) 
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1  THE  BGK  MODEL 


by  eliminating  time  derivatives  of  p,  u  and  T  in  favor  of  spatial  derivatives  using  the  zero 
order  (Euler)  form  of  Eqs.  (26-27).  The  final  expression  reads (Ferziger  and  Kaper,  1972): 


/(l)  =  _/(0)i  f 


5\  1  1  c2 

-J  coVr  log  T  +  —cc:  Vru-  Vro  U 


Inserting  expression  (35)  into  Eqs.  (28)  provides  Navier-Stokes  and  Fourier  law  closures 
for  stress  tensor  and  heat  flux,  respectively.  More  precisely  the  following  expressions  are 
obtained: 


P(1)  —  mj  ccf^dc  =  —2-S  (36) 

qW  =  Jj  R^VrT  (37) 

where  S  is  the  traceless  symmetric  component  of  Vu.  As  is  clear  from  the  equations  above, 
the  BGKW  model  predicts  the  following  expressions  for  the  shear  viscosity  coefficient  p 
and  thermal  conductivity  A: 

»=P-,  A  M  (38) 

v  2  v 

In  principle,  Eqs.  (38)  provide  an  easy  recipe  to  fit  the  single  model  adjustable  quantity, 
v(p,  T )  to  the  properties  of  a  specific  gas.  However,  it  can  be  readily  verified  that  the 
BGKW  model  predicts  a  unit  value  of  the  Prandtl  number  Pr  =  whereas  the  correct 
value  is  slightly  temperature  dependent  and  close  to  2/3  for  a  monatomic  gas.  Therefore, 
the  velocity  independent  collision  frequency  v  can  be  tuned  to  obtain  either  the  correct 
gas  viscosity  or  the  correct  gas  thermal  conductivity  but  not  both. 

Before  describing  and  discussing  the  model  modifications  and  extensions,  which  aim  at 
eliminating  drawbacks  and  extend  its  applications  to  more  complex  fluids,  it  is  worth 
considering  a  few  test  problem  where  the  flow  properties  obtained  by  the  model  described 
above  are  compared  with  the  solutions  of  the  full  Boltzmann  equation. 

1.3  BGKW  Model  Applications 

The  properties  of  the  BGKW  kinetic  model  with  velocity  independent  collision  frequency 
are  discussed  below  on  the  basis  of  the  following  simple  test  problems: 

1.  Normal  shock  profiles  in  a  hard  sphere  gas 

2.  1-D  heat  flow  in  a  rarefied  hard  sphere  gas  confined  between  parallel  plates. 

3.  Low  Mach  number  2D  flow  in  a  driven  square  cavity 

1.3.1  Normal  shock  profiles  in  a  hard  sphere  gas 

The  propagation  of  a  planar  shock  wave  is  studied  in  the  wave  front  reference  frame.  The 
resulting  stationary  flow  field  is  assumed  to  be  governed  by  the  one-dimensional  steady 
BGKW  equation 


~ f) 


(39) 
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1.3  BGKW  Model  Applications 


x  being  the  spatial  coordinate  which  spans  the  direction  normal  to  the  (planar)  wave 
front.  The  model  collision  frequency  is  computed  as 

=  Lifco’  '‘hs<t)  =  i^w^  (40) 

being  p  =  pRT  the  pressure  and  Phs(T)  the  first  approximation (Ferziger  and  Kaper, 
1972)  of  the  viscosity  of  a  gas  of  hard  spheres  of  mass  m  and  diameter  a.  It  is  further 
assumed  that,  far  from  the  wave  front,  the  distribution  function  f(x,  v )  satisfies  the 
boundary  conditions 


lim  f(x,  v)  =  <hT(u) 

x — ^-^poo 


rr 


(2ttRT^)3/2 


exp 


V*)2  +  v2y  +  v2z 
2  RT^ 


(41) 


where  uT,  VT  and  TT  are  the  upstream  and  downstream  values  of  number  density,  velocity 
and  temperature,  respectively.  The  parameters  of  the  equilibrium  states  specified  by  Eq. 
(41)  are  connected  by  the  Rankine-Hugoniot  relationships 

w  =  rp  =  MM-)2  T+  =  [5(A/~)2  —  1]  pf-)2  +  3] 

V+  n-  (M-)2  +  3  T~  16(M-)2  V  ’ 

In  Eqs.  (42)  M~  denotes  the  upstream  infinity  Mach  number  defined  as 


M~ 


V~ 

(-yRT-)'12 


(43) 


being  7  =  5/3  the  specihc  heat  ratio  of  a  monatomic  gas. 

Shock  profiles  obtained  from  the  numerical  solution  of  Eq.  (39)  are  compared  to  the  full 
Boltzmann  results,  obtained  from  DSMC  simulations,  in  Figures  1  and  2.  As  expected, 
the  faster  thermalization  intrinsic  in  the  BGKW  approximation  produces  steeper  density 
profiles.  The  advance  of  the  temperature  profile  rise  is  also  smaller,  although  a  longer 
upstream  tail  is  present  in  the  M~  =  5.0.  Moreover  the  slight  temperature  overshoot 
present  in  BE  solution,  seems  to  be  absent  in  the  BGKW  solution. 


1.3.2  1-D  heat  flow  in  a  rarefied  hard  sphere  gas 

The  limitations  on  the  accuracy  of  BGKW  model  predictions  imposed  by  the  incorrect 
value  of  the  Prandtl  number  clearly  appear  in  the  study  of  the  one-dimensional  flow  of  a 
monatomic  gas  confined  in  the  gap  between  two  infinite  parallel  plates  located  at  points 
x  =  ±L/2  of  a  reference  frame  in  which  the  x  coordinate  spans  the  direction  normal  to  the 
plates,  kept  at  constant  and  uniform  temperatures  Twl  and  Tw2 ,  respectively.  Without  loss 
of  generality,  it  is  assumed  that  the  ’’cold”  plate  position  is  x  =  —L/2  and  its  temperature 
is  Tw  1  whereas  the  ’’hot”  plate  position  is  x  =  +L/2  and  its  temperature  is  Tw2  >  Twl. 
The  gas  motion  is  assumed  to  be  governed  by  Eq.  (39)  in  which  the  collision  frequency  is 
obtained  by  Eqs.  (40).  The  boundary  conditions  at  the  wall  are  assigned  adopting  a  simple 
Maxwell’s  gas  surface  interaction  model  which  mixes  diffuse  and  specular  re-emission  with 
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Shock  Profiles  from  BE,  BGKW  and  ES  Model 

Hard  Sphere  Gas  -  M  ^=2.0 


xfkl 


Figure  1:  Normal  shock  profiles  in  a  hard  sphere 
gas,  M~  =  2.0.  Black:  normalized  density  n{x)/n~\ 
green:  normalized  temperature  T(x)/T~ ;  red:  nor¬ 
malized  mean  velocity  (x-component)  ux(x) / \J RT~ . 
Solid  lines:  DSMC;  symbols:  direct  solution  of  the 
Boltzmann  equation;  dashed  lines:  BGKW  model 


weights  cx  and  1  —  cx,  respectively.  Accordingly: 


f/  L  \  IT'wl 

n~rv)  =  aJ^mi^exp 
Hrv)  =  a(2,flr„)^exp 


V- 


2  RTwl 

.2 


v 


2  AT, 


w  2 


+  (!  W<0 

+  (!  vx>0 


-vx,vy,vz) 


v  = 

Mass  conservation  at  boundaries  determines  the  densities  of  wall  Maxwellians  as 


2tt 

RTwi 


'  VX<Q 


vxf(—^,v)dv 


2tt 


AT, 


w  2 


'  vx>0 


Vxf(^r,v)dv 


(44) 

(45) 

(46) 


(47) 


It  is  esily  found  that  the  solutions  of  the  problem  depend  on  the  three  parameters  ex, 
T„,2 / Tw i  and  Kn.  The  Knudsen  number  Kn  is  defined  as  Ao/T,  being  Ao  =  a  ref" 

erence  value  for  the  mean  free  path  obtained  from  the  mean  density  no  =  ^  S-l/2  n(x)  d x ■ 
Density,  temperature  and  heat  flux  profiles  obtained  by  numerical  solutions  of  Eq.  (39) 
are  compared  to  DSMC  simulations  of  the  hard  sphere  dilute  gas  in  Figure  3  in  the  case 
a  =  0.826,  Tw2/Twi  =  1.14,  Kn  =  0.7582.  It  is  worth  observing  gas  rarefaction  produces 
the  typical  temperature  jumps  at  walls  where  the  gas  temperature  differs  from  the  wall 
temperature.  The  agreement  between  the  the  DSMC  and  BGKW  density  and  tempera¬ 
ture  profiles  is  quite  good.  However,  as  expected,  a  discrepancy  is  found  on  the  heat  fluxes 
predicted  by  the  two  models  because  the  BGKW  collison  frequency  has  been  obtained 


RTO-EN-AVT-194 


3-9 


1  THE  BGK  MODEL 


1.4  Driven  Cavity  Flow 


Shock  Profiles  from  BE  and  BGKW  model 

Hard  Sphere  Gas  Mm=5.0 


Figure  2:  Normal  shock  profiles  in  a  hard  sphere 
gas,  M~  =  5.0.  Black:  normalized  density  n{x)/n~ ; 
green:  normalized  temperature  T(x)/T~ ;  red:  nor¬ 
malized  mean  velocity  (x-component)  ux(x)/\ / RT~ . 
Solid  lines:  DSMC;  dashed  lines:  BGKW  model 


from  viscosity.  A  more  complete  data  set  is  given  in  Table  ??  which  reports  values  of  the 
ratio  of  the  computed  heat  flux  qx  to  the  heat  flux  value  in  free  molecular  flow  conditions 
(Kn  — >  oo),  qxM ,  as  function  of  the  Knudsen  number  Kn  for  a  =  0.826,  Tw2/Tw\  =  1.14. 
As  is  clear,  tuning  collision  frequency  on  viscosity  leads  to  incorrect  predictions  of  heat 
fluxes  but  the  accuracy  can  be  recovered  by  obtaining  v  from  the  correct  thermal  con¬ 
ductivity  A hs  °f  the  hard  sphere  gas.  Of  course,  such  choice  would  spoil  any  shear  stress 
prediction. 

1.4  Driven  Cavity  Flow 

The  two  test  problems  described  above  show  that  the  simpler  mathematical  structure 
of  the  BGKW  model  is  not  without  consequences  on  the  accuracy  of  rarefied  flow  pre- 


Kn 

Boltzmann  Equation1 

BGKW  (u  from  /i)2 

BGKW  (u  from  A)1 

0.7582 

0.7558 

0.6759 

0.7566 

0.2994 

0.5807 

0.4774 

0.5756 

0.1942 

0.4843 

0.3811 

0.4786 

0.1395 

0.4094 

0.3117 

0.4041 

0.0658 

0.2538 

0.1813 

0.2512 

Table  1:  Heat  flux  ratio  qx/qxM  as  a  function  of  Kn;  a  =  0.826,  Tw2/Tw i  =  1.14 
;  1  -  (Ohwada,  1996);  2  -  (Graur  and  Polikarpov,  2009). 
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1,14 
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Figure  3:  Heat  flux  between  parallel  plates.  Hard  sphere  gas,  a  =  0.826,  TW2/TW i  =  1.14, 
Kn  =  0.7582. 


1  i  1  i  1  i  r 

-  qxU)/q0  (BE) 

--  qjx)/qn  (BGKW) 


_L 


_L 


_L 


_L 


_L 


_L 


_L 


-0,6  -0,5  -0,4  -0,3  -0,2  -0,1 


x/A,n 


0,1  0,2  0,3  0,4  0,5  0,6 


Heat  flux  between  parallel  plates 

Comparison  of  BE,  ES  and  BGKW  kinetic  models 


dictions.  However,  there  are  situations  where  the  application  of  typical  numerical  tools 
used  to  obtain  solutions  of  the  full  Boltzmann  equation  becomes  difficult.  For  instance, 
in  standard  implementations  of  DSMC  schemes  it  is  difficult  to  reduce  statistical  noise  to 
describe  accurately  the  tiny  deviations  from  equilibrium  met  in  low  Mach  number  micro¬ 
flows.  In  such  situations,  the  adoption  of  kinetic  models  allows  adopting  deterministic 
schemes  which  can  easily  capture  small  deviation  from  equilibrium.  The  following  exam¬ 
ple  shows  that  sensible  tuning  of  the  BGKW  model  collision  frequency  produces  a  fairly 
accurate  description  of  two-dimensional  rarefied  gas  flows. 

A  monatomic  gas  is  confined  in  the  two-dimensional  square  cavity 


C  =  {(x,y)  :  0  <  x  <  L,  0  <  y  <  L} 


The  flow  is  driven  by  a  uniform  translation  of  the  top  with  velocity  Vwex ,  being  ex  a  unit 
vector  parallel  to  x  direction.  The  gas  flow  is  governed  by  the  two-dimensional  steady 
BGKW  equation 

+ =  "(4>  “ f)  (48) 

where  the  collision  frequency  is  obtained  from  the  viscosity  of  the  hard  sphere  gas,  as 
discussed  above.  It  is  further  assumed  that  all  the  walls  are  kept  at  uniform  and  constant 
temperature  Tw  and  that  the  gas  atoms  which  strike  the  walls  are  re-emitted  according 
to  the  Maxwell’s  scattering  kernel  with  complete  accommodation 


f(x,v) 


nw  (s?) 

(2nRTwy/*  6XP 


[v  -  Vw(x)]2  1 
2  RTW  J’ 


( v  —  Vw)  o  h  >  0 


(49) 


where  a?  is  a  point  of  the  boundary,  h(x)  the  inward  normal  at  x,  Vw(x)  is  the  wall 
velocity,  different  from  zero  only  on  the  top  wall,  Tw  the  wall  temperature  and  nw(x)  the 
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wall  density  which  is  determined  by  impinging  mass  flux  through  the  following  relationship 

/  2n  \  1 /2  f 

nw(x)  =  -==-  /  |(v  -  Vw)  o  n\f  dv  (50) 

\JX1WJ  J(v—Vw)on<0 

which  ensures  zero  net  mass  flux  at  boundary  points. 

It  is  assumed  that  the  initial  gas  state,  at  time  t  —  0,  is  described  by  the  uniform  equilib¬ 
rium  Maxwellian 

/(x,  l>|0)  =  *„(«<)  =  (2^0)3/2  eXP  (" ^jr)  (51) 

being  no  and  T0  =  Tw  the  initial  values  of  the  uniform  density  and  temperature,  respec¬ 
tively.  The  non-dimensional  form  of  the  governing  equation  is  easily  obtained  by  adopting 
the  reference  mean  free  time  To  =  l/z^o  and  mean  free  path  Ao  =  y/2RT0ro  as  time  and 
length  units,  respectively.  It  is  also  immediately  seen  that  the  problem  solutions  depend 
on  the  dimensionless  wall  velocity  Vw/ y/2RTo  and  the  rarefaction  parameter  5  =  L/ Ao 
which  is  the  reciprocal  value  of  the  Knudsen  number  Kn  =  Ao / L. 

Eq.  (48)  can  be  easily  solved  by  a  variety  a  finite  difference  schemes(?).  The  BGKW 
velocity  field  is  compared  to  the  full  Boltzmann  equation  results  in  Figure  4.  It  should  be 
noted  that  the  top  wall  velocity  ratio  Vw/\]2 RT0  has  been  set  equal  to  10-2  and  that  the 
full  Boltzmann  equation  has  been  solved  by  a  semi-deterministic  method  which  adopts 
a  Monte  Carlo  quadrature  scheme  to  compute  the  Boltzmann  collision  integral  given  in 
Eq.  (2).  As  Figure  4  shows,  the  agreement  of  BGKW  and  BE  flow  fields  is  pretty  good. 
Moreover,  BGKW  computing  time  is  a  small  fraction  of  the  computing  time  of  Boltzmann 
equation  solutions.  Hence,  kinetic  models  are  of  some  interest  for  micro-flows  modeling. 
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Figure  4:  Profiles  of  the  dimensionless  (a)  horizontal  mean  velocity 
along  the  vertical  line  crossing  the  center  of  the  cavity,  and  (b)  ver¬ 
tical  mean  velocity  component  along  the  horizontal  line  crossing  the 
center  of  the  main  vortex.  Solid  and  dashed  lines:  numerical  solutions 
obtained  by  solving  the  full  Boltzmann  equation  with  a  semi-regular 
method (Frezzotti  et  al.,  2010)  for  8  —  10  and  8  =  0.1,  respectively.  Cir¬ 
cles  and  squares:  numerical  solutions  reported  in  Ref.  (Varoutis  et  ah, 
2008)  for  8  =  10  and  8  —  0.1,  respectively. 
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2  BGKW  model  extensions 

Several  attempts  have  been  made  to  improve  the  BGKW  model  and,  in  particular,  to 
obtain  kinetic  models  that,  while  keeping  its  simplicity,  do  have  the  correct  hydrodynamic 
limit.  As  far  as  the  linearized  form  of  the  kinetic  equation  is  concerned,  there  exist  a 
systematic  procedureCercignani  (1988)  to  construct  a  sequence  of  linear  models  which 
contains  the  linearized  BGKW  model  as  first  element  and  it  has  the  linearized  Boltzmann 
operator  as  final  limit.  For  some  linearized  models  a  non-linear  counterpart  can  be  found. 

2.1  The  Ellipsoidal  Statistical  Model 

In  1963  Holway(Holway,  1963)  and,  independently,  Cercignani  proposed  the  Ellipsoidal 
Statistical  model  (ES-model)  which  keeps  the  overall  structure  of  the  BGKW  model,  but 
approximates  the  gain  term  in  Eq.  (1)  as  the  product  of  the  collision  frequency  v  and  an 
anisotropic  Gaussian.  More  precisely,  the  ES-model  kinetic  equation  has  the  form: 

^  +  v-Vrf  =  u[G(r,v\t)  —  f(r,v\t)}  (52) 

G(r,v\t)  =  >l^  =  exp  1  :  cc j  (53) 

y(27r)3detA  '  ' 

In  Eq.  (57)  c  =  v  —  u(r\t)  is  the  peculiar  velocity,  whereas  the  tensor  A  is  defined  as: 

A  =  (l-  (3)RTI  +  /3-P  (54) 

P 

being  /3  =  — l—-  It  can  be  shown  that  the  collision  term  of  the  ES-model  conserves 
mass,  momentum  and  energy  and  leads  to  a  hydrodynamic  limit  with  a  specified  Prandtl 
number  equal  to  Pr,  provided  the  relationship  between  collision  frequency  and  viscosity 
is  modified  as: 

v  =  Pr—  (55) 

P 

Unfortunately,  ES-model  more  complicated  structure  makes  the  study  of  some  of  its  prop¬ 
erties  more  difficult.  For  instance,  the  proof  of  the  77-theorem  for  ES-model  is  far  from 
being  trivial  and  it  has  been  given  only  a  few  years  ago(Andries  et  ah,  2000).  On  the  other 
hand,  numerical  solutions  of  Eq.  (57)  can  be  easily  obtained  by  the  same  deterministic 
schemes  developed  for  the  BGKW  model  or  by  particle  schemes (Andries  et  ah,  2002b). 
Although  the  model  provides  the  correct  hydrodynamic  limit,  it  is  not  guaranteed  that 
it  will  yield  accurate  results  in  the  transition  regime.  Therefore,  an  investigation  of  the 
model  properties  through  the  study  of  simple  test  problems  is  advisable.  Normal  shock 
profiles  can  be  easily  obtained  by  using  Eq.  (57)  in  place  of  Eq.  (4)Ferziger  and  Kaper 
(1972).  Figure  5  compares  BGKW,  ES-model  and  full  Boltzmann  density  and  tempera¬ 
ture  profiles  in  a  stationary  shock  wave  in  a  hard  sphere  gas  with  upstream  infinity  Mach 
number  =  2.0.  It  is  interesting  to  observe  that  the  slope  of  the  density  profile  in 
the  central  region  of  the  shock  is  closer  to  the  Boltzmann  equation  results.  The  overall 
agreement  with  the  temperature  profile  is  also  generally  better  than  the  BGKW  model, 
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2.2  Shakhov  Model 


Shock  Profiles  from  BE,  BGKW  and  ES  Model 

Hard  Sphere  Gas  -  M  m=2.0 


Figure  5:  Normal  shock  profiles  in  a  hard  sphere  gas, 
M~  =  2.0. 


however  the  ES-model  exhibits  a  more  pronounced  temperature  tail  which  extends  up¬ 
stream. 

As  shown  in  Table  2.1,  the  improvements  brought  by  the  ES-model  are  also  evident  in  the 
heat  flux  problem  described  in  the  previous  section  and  recently  studied  in  Ref.(Graur 
and  Polikarpov,  2009)  on  ground  of  different  kinetic  models.  Es-model  applications  to 
hypersonic  flows  have  been  described  in  Ref.Andries  et  ah  (2002b)  which  also  contains 
the  description  of  a  particle  scheme  for  its  numerical  solution. 


2.2  Shakhov  Model 

As  is  clear  from  the  expression  of  the  gain  term  in  Eq.  (57),  the  strategy  to  extend  BGKW 
model  capabilities  is  to  replace  a  Maxwellian  with  function  containing  a  larger  number  of 
adjustable  parameters.  The  Gaussian  function  is  by  no  means  the  unique  possible  choice. 
In  order  to  obtain  a  model  with  the  correct  Prandtl  number,  Shakhov  (Shakhov,  1968) 


Kn 

QJ€M  (a) 

QM'M  (b) 

0.7582 

0.7558 

0.7488 

0.2994 

0.5807 

0.5684 

0.1942 

0.4843 

0.4719 

0.1395 

0.4094 

0.3980 

0.0658 

0.2538 

0.2466 

Table  2:  Heat  flux  ratio  qx/qxM  as  a  function  of  Krr,  a  =  0.826,  Tw2 /Tw\  =  1.14 
;  (a)  Boltzmann  Equation(Ohwada,  1996);  (b)  ES-model  ( u  from  n,  Pr  =  2/3)(Graur 

and  Polikarpov,  2009) 
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proposed  the  following  model  (S-model): 


y^  +  v-X7rf  =  u[S(r,v\t)~  f(r,v\ 

Pr 


S(r,v\t )  =  $(ryu|f) 


1  + 


5  pRT 


co  q 


(56) 

(57) 


S-model  gain  term  is  the  product  of  the  local  Maxwellian  times  a  third  degree  polynomial 
of  peculiar  velocity  components.  Although  it  generally  performs  better  than  BGKW 
medcl,  it  should  be  observed  that  a  proof  of  the  (iL)-theorem  is  still  lacking  for  S-model. 
Moreover,  the  polynomial  which  multiplies  may  become  negative  in  some  regions  of 
the  velocity  space.  Since  negative  values  usually  occur  in  high  velocity  regions,  negative 
values  of  the  gain  term  do  not  harm  in  modeling  low  speed  flows  where  the  model  has 
been  successfully  applied(Graur  and  Polikarpov,  2009). 


2.3  Models  with  velocity  dependent  v 

Kinetic  models  with  the  correct  Prandtl  number  can  also  be  obtained  by  allowing  the  col¬ 
lision  frequency  v  to  depend  on  molecular  velocity  v ,  as  it  does  in  the  Boltzmann  equation 
where  the  collision  frequency  is  a  functional  of  /,  which  for  hard  sphere  interaction,  takes 
the  form  shown  in  Eq.  (2) 


v[r,v)  =  °—  /  dv i  /  d2k  f(x,Vi\t)\k  ■  vr\  =  n(r)ira2vr  (58) 

2  Jr3  Js 2 

As  shown  in  Ref.  (Mieussens  and  Struchtrup,  2004),  when  the  collision  frequency  v  in 
Eq.  (4)  depends  on  velocity,  then  the  Prandtl  number  is  given  by  the  following  expression: 


Pr 


2R\ 


I  exp  (-Tj2)  dr] 

f  exP  (~d2)  dV 


(59) 


where  r/  =  c/\/2RT.  The  model  greater  generality  is  obtained  at  the  cost  of  some  com¬ 
plications.  Actually,  the  macroscopic  fields  n(r\t),  ii(r\t),  T(r\t)  no  longer  coincide  with 
local  density,  bulk  velocity  and  temperature.  The  comprehensive  study  in  Ref.  (Mieussens 
and  Struchtrup,  2004)  also  show  that  frequency  dependent  model  do  allow  some  improve¬ 
ments  if  a  sensible  choice  for  u(c)  is  made.  However,  the  choice  is  problem  dependent  as 
axpected  on  ground  of  expression  (58) 
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3  Kinetic  models  for  mixtures 


The  aim  of  this  section  is  to  briefly  sketch  the  methods  and  models  which  lead  to  the 
extension  of  the  kinetic  models  described  above  to  gaseous  mixtures.  For  simplicity, 
the  case  of  a  binary  mixture  of  monatomic  species  in  absence  of  external  fields  will  be 
considered. 

The  mixture  is  described  by  a  system  of  two  coupled  Boltzmann  equation  in  the  form: 


2 


df 

3= 1 


The  collision  integrals  Dq  (/*,/))  satisfy  the  following  relationships: 


j  CM,  fj)A(v)  dv  =  0,  ^(v)  =mi,miV,-miV2,i  =  1,2 

/c, //,,/,)</«  =  o,  i  /  j 

m  I  fj)vdv  +  mj  f  fj)v  dv  =  0,  i  r  / 

">i  f Cijdi,  fj)v2dv  +  mj  f  Cjilfi,  fj)v2  dv  =  0,  i  /  / 


(60) 


(61) 

(62) 

(63) 

(64) 


expressing  conservation  of  mass  momentum  and  energy  in  self- interaction  of  individual 
species,  conservation  of  mass  in  cross-interactions  and  conservation  of  total  momentum 
and  energy  in  cross-interaction. 

Following  the  arguments  leading  to  the  BGKW  model,  it  is  quite  natural  to  replace  each 
collision  term  in  the  Boltzmann  equations  (60)  with  a  relaxation  term: 


dfi 

dt 


+  V  ■  Wrfi 


2 

Uij  ij  ~  fi) 

3= 1 


nij(r\t) 

- - - exp 

[InR^irlt)}312 


[v  -  Uij{r\t)f  1 
2RiTij(r\t)  J 


(65) 

(66) 


It  is  immediately  realized  that  the  twenty  four  model  disposable  parameters/fields  vl3 , 
nVJ ,  TtJ  and  ui3  cannot  be  determined  by  the  the  sixteen  equations  61  and  the  request 
that  the  total  number  of  cross-collisions  is  (obviously)  the  same  for  both  components: 

n.iUji  =  njVij  (67) 

There  is  no  unique  way  to  define  the  quantities  ,  n%3 ,  T%3  and  u%3.  Hamcl(Hamel, 
1965),  for  instance,  adds  to  Eqs.  (61,67)  the  request  that  the  expressions  for  momen¬ 
tum  and  energy  transfers  in  cross-collisions  is  the  same  as  in  a  mixture  of  Maxwell’s 
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molecules (Cercignani,  1988).  The  resulting  expressions  read(Hamel,  1965): 


Uu  ( 

>|f) 

=  ni:j(r\t )  =  ni(r\t )  = 

j  fi(r,v\t )  dv 

(68) 

Uii\ 

>1 1) 

=  Ui{r\t)  =  -  [  vfi(r, 
n  J 

v\t)  dv 

(69) 

Uij  ( 

>1 1) 

=  Ui(r\t)  +  aij(ui(r\t) 

~ Uj(r\t )) 

(70) 

Ta\ 

tI  t) 

=  Ti(r|t)  =  3  urJ(V 

-  Ui)2fi(r,  v\t)  dv 

(71) 

Tu( 

>1 1) 

=  r,(r|f)  +  fiij[Tj(r\t)  - 

-  Ti(r\t )]  +  7 ij[ui(r\t)  -  Uj(r\t )]2 

(72) 

being  al3 ,  and  7 ^  constants  depending  on  molecular  masses  and  interaction  potentials 
strenght.  In  the  case  of  Maxwell’s  molecules  collision  frequencies  can  be  written  as  vl3  = 
rijKij ,  being  k13  interaction  potential  constants  that  can  be  determined  to  match  individual 
components  an  mixture  viscosities. 

It  is  to  be  observed  that  a  more  systematic  procedure  to  obtain  kinetic  models  for  mixtures 
is  possible  for  linearized  kinetic  models(Mc  Cormack,  1973). 

It  should  also  be  observed  that  BGKW-like  kinetic  model  for  mixtures  in  the  form  of 
Eqs.  (65)  generally  do  not  reduce  to  one-component  model  when  one  considers  a  fictitious 
mixture  of  mechanically  identical  components  (Garzo  et  al.,  1989).  Consistent  kinetic 
models  which  correct  such  inconsistency  have  been  proposed  in  Refs. (Garzo  et  ah,  1989; 
Andries  et  ah,  2002a). 
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4  Kinetic  models  for  polyatomic  gases 

The  problem  of  bringing  the  effects  of  internal  molecular  structure  into  kinetic  equations 
is  not  trivial  for  several  reasons.  Hence,  the  present  notes  does  not  even  attempt  to  give  a 
complete  account  of  kinetic  theory  of  dilute  polyatomic  gases  which  is  properly  described 
in  a  few  textbooks.  It  will  simply  be  noted  that  the  difficulty  of  describing  collisional 
processes  between  molecules  with  internal  degrees  of  freedom  (rotational,  vibrational, 
electronic)  reduces  the  ’’distance”  between  a  description  based  on  the  full  Boltzmann 
equation  and  kinetic  model.  As  a  matter  of  fact,  the  former  is  sometimes  based  on  quite 
phenomenological  collision  models.  The  content  of  the  following  notes  is  therefore  limited 
to  the  short  presentation  of  the  widely  used  Bognakke- Larsen  collision  model  and  of  one 
of  the  possible  extensions  of  the  BGKW  kinetic  model.  The  predictions  of  the  two  poly¬ 
atomic  gas  models  will  be  compared  on  a  simple  one-dimensional  flow  resulting  from  the 
intense  evaporation  from  a  planar  surface. 

The  problem  treatment  presented  here  is  limited  to  the  temperature  range  where  vibra¬ 
tional  excitation  can  be  neglected.  Moreover,  since  the  gap  between  quantized  rotational 
energy  levels  is  supposed  to  be  much  smaller  than  nTref  ( Tref  is  a  reference  characteristic 
flow  temperature  value),  vapor  molecules  are  supposed  to  behave  as  classical  rigid  rota¬ 
tors  of  average  diameter  a  and  mass  m  with  j  —  2  (linear  molecule)  or  j  =  3  (non-linear 
molecule)  rotational  degrees  of  freedom.  It  is  assumed  that  the  gas  motion  is  governed  by 
the  following  one-dimensional  Boltzmann  equationCercignani  (1988);  Kuscer  (1989): 

Zx^(v,£\x)  =  I  [f(v,1,£[\x)f(v,,£'\x)~f(v1,£1\x)f(v,£\x)}Q£^d3^1d£1  (73) 

being  f(v,£\x)  the  distribution  function  of  molecular  velocity  v  and  rotational  energy  S 
at  location  x.  It  is  to  be  noted  that  rotational  degrees  of  freedom  are  taken  into  account 
only  through  the  rotational  energy  £,  not  through  the  whole  set  of  angular  coordinates  and 
momenta.  In  Eq.(73),  Q  is  defined  as 

p  pE-S'  pE~£[  C'2 

Q  —  d2e!  /  £,fld£'  £g  d£[^a(E-  e!  o  e;£',  £[  -»•  £,  Si)  (74) 

Js  J  0  Jo  Sr 

being  a(E;e'  o  e;  £',  £[  — *  £,£f)  the  differential  cross-section  associated  with  a  binary 
collision  which  produces  a  pair  of  molecules  in  the  final  states  (v,£),  (v±,£i)  from  a  pair 
of  molecules  in  the  initial  states  (v',£'),  (v[,£[).  The  argument  E  denotes  the  conserved 
total  energy  in  the  center  of  mass  reference  frame: 

E  =  +  £  +  £i  =  \mg2  +  £'  +  £[  (75) 

The  unit  vectors  e  —  ff-  and  e  =  ^  have  the  directions  of  the  relative  velocities  v'  = 

Sr  Sr 

v\  —  v'  and  vr  =  tq  —  v  before  and  after  a  collision,  respectively.  The  exponent  //  in 
Eq.(73)  takes  the  values  0  for  j  =  2  and  |  for  j  =  3. 

4.1  The  collision  model 

The  collision  dynamics  and  cross-section  have  been  obtained  from  a  phenomenological 
model  proposed  by  Borgnakke  and  LarsenBorgnakke  and  Larsen  (1975).  The  model 
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overcomes  some  of  the  restrictions  and  complications  of  previous  mechanical  models  of 
translational-rotational  coupling  (rough  spheres,  loaded  spheres,  spherocylinders) Chapman 
and  Cowling  (1990)  since  it  can  be  easily  adapted  to  reproduce  experimental  translational- 
rotational  relaxation  rates  with  good  accuracy Wysong  and  Wadsworth  (1998).  Moreover, 
the  collision  algorithm  derived  from  the  model  is  very  well  suited  to  particle  schemes  used 
to  obtain  numerical  solutions  of  the  Boltzmann  equationBird  (1994).  In  the  particular 
form  of  the  Borgnakke- Larsen  model  adopted  here,  collision  dynamics  is  organized  as 
follows: 


•  The  collision  probability  of  two  molecules  in  the  pre-collision  states  (v',£'), 

is  proportional  to  where  ahs  =  n a2  is  the  integral  cross-section  of  hard  sphere 

molecules  and  £'  =  —  i/||  is  the  relative  velocity  modulus. 

•  An  individual  collision  is  inelastic  with  probability  z  or  clastic  with  probability 
1  —  z.  An  inelastic  collision  gives  rise  to  an  exchange  between  translational  and 
rotational  energies,  as  explained  below.  In  an  clastic  collision  pre-  and  post-collision 
rotational  energies  do  not  change,  i.e.  £  =  £',  £\  =  £[.  Conservation  of  total 
energy  then  implies  =  £'  and,  according  to  hard  sphere  impact  dynamics,  post¬ 
collision  relative  velocity  is  written  as  vr  =  £re,  being  e  a  random  vector  uniformly 
distributed  on  the  unit  sphere  S. 


•  In  an  inelastic  collision  total  energy  E  is  randomly  partitioned  between  translational 
and  rotational  motion  by  sampling  the  translational  energy  fraction  Etr/E  from 
a  given  probability  density  function  V\  (Etr / E\j) .  The  available  total  rotational 
energy  Erot  —  £  +  £i  —  E  —  Etr  is  then  randomly  distributed  between  the  collision 
partners  by  sampling  the  fraction  £ / Erot  from  a  given  probability  density  function 
V‘2(£ / Erot\j).  The  relative  velocity  after  a  collision  is  again  written  as  vr  =  £ re , 

where  e  is  a  random  unit  vector  and  £r  =  .  j 


The  specific  form  of  the  probability  densities  V\(£ / Erot\j)  and  V2(Etr/ E\j)  depends  both 
on  the  number  of  internal  degrees  of  freedom  and  on  the  assumed  intermolecular  interac- 
tionBird  (1994).  In  the  case  of  hard  sphere  interaction  and  j  —  2  they  take  a  particularly 
simple  formKuscer  (1989);  Bird  (1994) 


Vx{Etr!E\2) 
'E>2(£  /  Erot\2) 


(76) 

(77) 


As  shown  by  Eq.(76),  post-collision  translational  energy  has  a  parabolic  distribution;  the 
available  Erot  amount  is  then  randomly  divided  between  £  and  £i,  according  to  Eq.(77). 
Taking  into  account  the  assumed  scattering  isotropy  and  Eqs. (76,77),  the  collision  cross- 
section  takes  the  form: 

a(B;e'oe-,£',£;^£,£I)  =  (78) 

0(e,  ei,  e',  e))  =  (1  —  z)6(e  —  e')5(ei  —  e^)  +  6^(1  —  e  —  ei)  (79) 


being  e  =  £/E. 
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The  strength  of  translational-rotational  coupling  is  determined  by  the  mixing  param¬ 
eter  z  which  can  be  made  to  depend  on  the  local  flowficld  temperature  to  fit  experimental 
relaxation  ratesWysong  and  Wadsworth  (1998)  A  kinetic  model  which  closely  patterns 
BGKW  model  has  been  proposed  by  Holway  in  1966Holway  (1966).  As  usual,  relaxation 
terms  replace  the  full  Boltzmann  collision  integral: 

%  +  ^  =  "ei(n,  Tt)  ($e,  -  /)  +  uin(n,  Tt)  ($in  -  /)  (80) 

where  the  “frozen”  local  equilibrium  and  local  equilibrium  Maxwellian  distribution  func¬ 
tions,  $e/('U,£)  and  <3 ?in(v,£),  are  defined  as  follows: 


i  £) 


n(£\x,  t)  \ 

v  —  u(x ,  t)] 

[2tt  BTt(x,t)]W  1  \ 

2  RTt(x,t) 

N(x,t)  f 

[v  —  u(x,t )] 

[2tt RT(x,  f)]3/2  1  \ 

2 RT(x,  t ) 

gj/2-l 

£ 

r(j/2)[KT(x,i)p/2<ixp 

L  nT(x,t) 

(81) 


(82) 


In  Eq.(80)  a  BGK-likeBhatnagar  et  al.  (1954)  expression  replaces  the  complicated  collision 
integral  of  Eq.(73).  The  first  term  describes  clastic  collisions,  whereas  the  second  one 
models  inelastic  collisions.  The  macroscopic  fields  associated  with  <l>e;(v,£)  and  £) 

are  defined  as: 


n(£\x,t)  = 

J  f(v,£\x,t)  d3£ 

(83) 

N(x,t )  = 

j  f(v,  £ \x,  t)  d3£  d£  =  Jn(£\x,t)d£ 

(84) 

u(x,t)  = 

N^x  ^  J  vf(v,£\x,t)d?£d£  =  ux(x,t)x  +  uy(x,t)y 

(85) 

Tt(x,t )  = 

s R.N(x,t)S[v  u{xM 2^v-£'x’t:>'^i£ 

(86) 

Tr{x,t)  = 

.  J,  [  £f(v,£\x,t)d3(d£ 

(87) 

jnN{x,t )  J 

T(x,t)  = 

3  Tt(x,t)  +jTr(x,t) 

3+? 

(88) 

The  quantity  n(£\x,  t)  is  the  number  density  of  molecules  having  internal  energy  £,  N(x,  t ) 
is  the  total  number  density,  u(x,t )  is  the  mean  velocity  of  the  gas,  whereas  Tt(x,t), 
Tr(x,t)  and  T(x,t)  are  the  translational,  internal  and  overall  temperature,  respectively. 
The  clastic  and  inelastic  collision  frequencies,  uei(n,  Tt)  and  um(n,  Tt)  have  been  computed 
from  the  shear  viscosity  //  of  the  hard  sphere  gas  as: 


vd{N,Tt) 

=  (1  -  z)vtot(N,Tt) 

(89) 

vin(N,Tt) 

=  zutot(N,  Tt) 

(90) 

vtot{ N,  Tt) 

NnTt 

(91) 

AT) 
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4.2  A  simple  application 


The  numerical  solution  of  Eq.(80)  can  be  greatly  simplified  by  a  transformation  first 
considered  by  ChuChu  (1965).  Let  us  define  the  reduced  distribution  functions  F(£x\x,t), 
G(^x\x,t),  H(£x\x,t)  and  K(£x\x,t)  as: 


/  F(^x\x,t)\ 
G{^x\x,t) 
H{kx\x,t) 
\K(kx\x,t)J 


(  1  \ 

(e2 + a 

V  £  ) 


f(v,£\x,t)d£yd£zd£ 


(92) 


It  can  be  easily  verified  that  the  reduced  distribution  functions  defined  above  obey  the 
following  system  of  kinetic  equations: 


where 


(F\ 

(F\ 

(  fc1-f\ 

(  Fin~F\ 

d 

G 

,  t  d 

G 

Gel  -  G 

Gin  -  G 

dt 

H\ 

+  £x7T“ 

ox 

H 

= 

Hel-H 

+  Dn 

Hin-H 

W 

W 

\Kel-K) 

\Kin  -  K j 

(  Fei(£,x\x,t)\ 


( 


G ei  (£# 

Hei  (£#  |x,  t) 

i 

|«£j  t)  J 

\ 

(  Fin(£,X\x,t)^ 

/ 

Ginijix  |«£? 

V 

Uy(x,t) 


\ 


jftTr 

2 


1  \ 

Uy(x,t ) 

Uy(x ,  t)  +  2 RT(x,  t ) 

jnT(x,t)  7 

2  / 


N(x,  t ) 

\f2iiRTt 


exp 


N(x,  t ) 
V2ttRT 


exp 


(■Cx  Ux )  1 

2  AT/  J 


(6  ~  ~»x)2  1 

2AT  J 


(93) 


(94) 


(95) 


4.2  A  simple  application 


We  consider  the  steady  one-dimensional  flow  of  a  polyatomic  vapor  evaporating  from  an 
infinite  planar  surface  (interface)  kept  at  constant  and  uniform  temperature  Tw .  The 
surface,  located  at  x  —  0,  separates  the  condensed  phase,  which  occupies  the  half-space 
x  <  0,  from  the  vapor  phase  flowing  in  the  half-space  x  >  0.  The  coordinate  x  spans  the 
direction  normal  to  the  surface  in  a  reference  frame  at  rest  with  respect  to  the  interface. 
The  gas  motion  has  been  computed  bu  solving  both  Eq.  (73)  and  Eq.  (80).  Molecular 
exchange  processes  between  the  vapor  and  condensed  phase  across  the  interface  have  been 
described  by  the  following  boundary  condition  for  /  at  x  —  0: 


vxf(v,e\0,t)  =  vxfe(v,e)  + 


K(v,e\v1,e1)\vlx\f(v1,e1\0,t)dv1de1,  vx  >  0  (96) 


'vlx<0 


The  distribution  function  of  evaporating  molecules,  fe.  is  assumed  to  be  a  half-range 
Maxwellian: 


fe  (v,e)  =  <7e 


N„ 


(27 tRTu 


\  3/2 


exp 


v 


2  AT, 


J/ 2~1 


X 


W/2)(kTw)R2 


exp 


kT„ 


Vr  >  0 


(97) 
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The  velocities  of  vapor  molecules  re-emitted  into  the  gas  phase  after  interacting  with  the 
interface  are  obtained  by  the  following  scattering  kernel  K  which  describes  a  pure  diffusive 
re-emissionCercignani  (1988) : 


K(v,e\vi,ei) 


(1-^2y]kvexp(-^) 

eJ/2-1  (  e  \ 

xv{j/2)^Twyn exp  (,-^J 


(98) 


In  Eqs.  (97,98)  Nw  denotes  the  saturated  vapor  density  at  temperature  Tw,  ae  is  the 
evaporation  coefficient,  k  is  the  Boltzmann  constant,  R  denotes  the  ratio  n/m  and  T  is 
the  complete  Gamma  function.  It  is  further  assumed  that  the  far  downstream  vapor  state 
is  described  by  an  equilibrium  Maxwellian  distribution  function: 


foo{v,e) 


iVoo 

wm^exp 


(v  -  Upp)2 
2RT00 


ej/2~l 

xf mw^exv 


(99) 


As  shown  in  Figures  6  and  7  the  full  Boltzmann  equation  corresponding  to  Borgnakke- 
Larsen  model  and  Holway’s  kinetic  model  produce  very  close  results. 
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x/X 

W 


Figure  6:  Density,  velocity  and  mass  flux  profiles  in  evaporation  flow 
of  a  polyatomic  gas  with  j  =  3  rotational  degrees  of  freedom;  z=0.3 
Mqo  —  0.534.  Numerical  data  from  solution  of  Holway’s  kinetic  model 
are  represented  by  solid  lines.  DSMC  results  represented  by  symbols: 
o,  N(x)/Nw;  >,  ux(x)/y/RTw;  x,  J[x)  =  N (x)ux(x) / Nwy/ RTW 
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W 


Figure  7:  Temperature  profiles  in  evaporation  flow  of  a  polyatomic  gas 
with  j  =  3  rotational  degrees  of  freedom;  z=0.3,  M, ^  =  0.534.  Nu¬ 
merical  data  from  solution  of  Holway’s  kinetic  model  are  represented 
by  solid  lines.  DSMC  results  represented  by  symbols:  o,  Tt(x)/Tw;  >, 
T\\  (x)/Tw;  A,  T±(x)/TW;D,  Tr(x)/Tw;  x,  T(x)/Tw. 
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